Formation au carroyage et lissage spatial sur R

Kim ANTUNEZ - Julien PRAMIL

octobre 2022

Plan :

1. Introduction
2. Configurations
3. Lissage
4. Indicateurs sur zonage à façon

Introduction

Objectifs du TP

• En 2018, création du package btb (PSAR analyse urbaine, Arlindo Dos Santo et François Sémécurbe).
• Sa principale fonction, kernelSmoothing, permet de réaliser très facilement un carroyage et un lissage sur des données géolocalisées avec R.

À partir de données ponctuelles, nous allons apprendre en utilisant le langage R :

• À carroyer les informations.
• À réaliser des lissages :
  • de densité,
  • de moyennes,
  • des lissages de taux,
  • et des lissages quantiles.
• À calculer un indicateur sur une zone à façon (données ponctuelles et de carroyées).

Liens utiles

• Code de la formation
• Site web des supports de formation

Avertissements

Secret statistique

Avant toute diffusion auprès des partenaires, il faut bien veiller à respecter :

• le secret
  • primaire
  • secondaire
  • fiscal
• les conventions établies avec les fournisseurs des données

Qualité des cartes

Pour simplifier : on prend des libertés avec la sémiologie cartographique

Auteur : Timothée Giraud, auteur de la librairie mapsf

Système de projection

Nom	Description	Code EPSG
Lambert93	Système de projection officiel pour la métropole	2154
LAEA	Système de projection européen	3035
WGS84	GPS (utile pour utiliser Leaflet)	4326

PI, voici les systèmes de projection que vous pouvez régulièrement rencontrer pour la métropole :

Configurations

Chargement des librairies

Cinq librairies sont nécessaires pour ce TP.

• sf pour manipuler des fichiers spatiaux (importer des .shp, transformer des projections, et réaliser des géotraitements) ;

• mapsf pour réaliser des cartes dans RStudio ;

• mapview (reposant sur leaflet) pour réaliser des cartes interactives (fond de carte OpenStreetMap);

• btb pour le carroyage et lissage ;

• dplyr pour le traitement des données, en particulier l’agrégation géographique.

Remarque : Le choix de dplyr plutôt que data.table se justifie ici du fait de sa forte compatibilité avec les objets géomatiques.

Charger les librairies nécessaires

## Liste des librairies utilisées
packages <-  c("dplyr","sf","btb","mapsf","leaflet","mapview")

## Vérifier si la librairie est installée, si non l'installer, puis la charger
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE, quiet = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

Sur le SSPCloud :

• Utiliser le package aws.s3 pour charger les données stockées dans Minio ;
• Plus rapide qu’un chargement classique.

if (!require("aws.s3", character.only = TRUE)) {
      install.packages("aws.s3",repos = "https://cloud.R-project.org",
                       dependencies = TRUE, quiet = TRUE)
      library("aws.s3", character.only = TRUE)
    }

Chargement de la base

• On utilise le base « Demandes de Valeurs Foncières »,
• Produite par la Direction générale des finances publiques (actes notariés).

Elle recense :

• l’ensemble des ventes de biens fonciers,
• Au cours des 5 dernières années,
• Hors Mayotte et en Alsace-Moselle.
• Biens bâtis ou terrains

À partir de cette source, nous avons constitué une base :

• Uniquement sur le périmètre de la petite couronne parisienne
• Pour l’année 2021.

Avec les 8 variables suivantes :

• id_mutation : identifiant unique de la vente
• date_mutation : date de la vente
• type_local : appartement ou maison
• nombre_pieces_principales : nombre de pièces dans le logement
• valeur_fonciere : prix de vente
• surface_reelle_bati : surface du logement
• x : longitude (en projection Lambert 93)
• y : latitude (en projection Lambert 93)

Remarques :

• Code de production de la base.

• Base non représentative de la réalité : usage purement pédagogique.

Dans le SSPcloud :

Importation de la base ventesImmo_couronneParis.RDS, stockée sous Minio

dans le “bucket public” : s3/projet-formation/r-lissage-spatial/.

# Charger la source de données (variable nommée dfBase) avec S3
url_bucket <- "https://minio.lab.sspcloud.fr/"
bucket <- "projet-formation"
object <- "r-lissage-spatial/ventesImmo_couronneParis.RDS"

dfBase <-
  aws.s3::s3read_using(
    FUN = base::readRDS,
    object = object,
    bucket = bucket
    ,
    opts = list("region" = "")
  )

En dehors du SSPCloud : téléchargement via URL.

url_file <- url(paste0(url_bucket,bucket,"/",object))
dfBase <- readRDS(url_file)

On peut ensuite manipuler notre base de données chargée en mémoire.

head(dfBase)

  id_mutation date_mutation  type_local nombre_pieces_principales
1 2021-447023    2021-01-08 Appartement                         3
2 2021-447024    2021-01-05 Appartement                         2
3 2021-447025    2021-01-08 Appartement                         3
4 2021-447027    2021-01-08 Appartement                         3
5 2021-447029    2021-01-05 Appartement                         2
6 2021-447030    2021-01-15 Appartement                         3
  valeur_fonciere surface_reelle_bati        x       y
1          480000                  64 647357.3 6868635
2          345000                  43 644483.5 6867695
3          384000                  41 648001.8 6866153
4          261900                  44 647414.8 6868937
5          407200                  24 646929.9 6864730
6         1040000                  90 645132.7 6864646

dim(dfBase)

[1] 34489     8

Chargement et préparation des données cartographiques

Fond de carte du territoire à étudier

• Chargement contour géographique des départements de la petite couronne parisienne.
• => Couches dites « vectorielles »
• Type de fichiers : .gpkg ou .shp
• On recommande le .gpkg (voir détails ici)

En cartographie, ces fichiers sont appelés des « couches vectorielles » (à ne pas confondre avec les couches « raster » qui correspondent, par exemple, aux fonds de carte satellite OpenStreetMap ou Google Maps). - Un .gpkg est un fichier unique alors qu’un fichier .shp tient dans 5 fichiers séparés et interdépendants. - Le format Géopackage est libre et ouvert, contrairement à Shapefile. - Un Shapefile impose des noms de variables de moins de 10 caractères, ce qui peut être pénible en pratique. - Le géopackage est devenu le format standard dans Qgis depuis la version 3 (logiciel libre de cartographie préconisé à l’Insee).

Récupérons donc le contour du territoire étudié avec la fonction st_read du package sf.

chemin_file <- paste0(url_bucket,bucket,"/r-lissage-spatial/depCouronne.gpkg")
depCouronne_sf <- sf::st_read(chemin_file)

Reading layer `depCouronne' from data source 
  `https://minio.lab.sspcloud.fr/projet-formation/r-lissage-spatial/depCouronne.gpkg' 
  using driver `GPKG'
Simple feature collection with 4 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
Projected CRS: RGF93 / Lambert-93

Visualisons cette couche vectorielle

head(depCouronne_sf)

Simple feature collection with 4 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 637307.1 ymin: 6843303 xmax: 671686 ymax: 6879253
Projected CRS: RGF93 / Lambert-93
  code           libelle reg surf                           geom
1   75             Paris  11  105 MULTIPOLYGON (((660897 6860...
2   92    Hauts-de-Seine  11  176 MULTIPOLYGON (((648796 6847...
3   93 Seine-Saint-Denis  11  237 MULTIPOLYGON (((659428 6861...
4   94      Val-de-Marne  11  245 MULTIPOLYGON (((656908 6846...

plot(depCouronne_sf$geom)

On renomme la variable geom en geometry

depCouronne_sf <- depCouronne_sf %>% rename(geometry=geom)

Une fois les départements chargés, nous allons sélectionner le contour de la commune de Paris.

# Sélection de Paris
paris_sf <- depCouronne_sf[depCouronne_sf$code=="75",]

Visualisons cette couche vectorielle

head(paris_sf)

Simple feature collection with 1 feature and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 643076 ymin: 6857499 xmax: 660897 ymax: 6867034
Projected CRS: RGF93 / Lambert-93
  code libelle reg surf                       geometry
1   75   Paris  11  105 MULTIPOLYGON (((660897 6860...

plot(paris_sf$geometry)

On reprojète dans le même système de projection (si besoin). Ici, en Lambert 93 (epsg 2154)

depCouronne_sf <- sf::st_transform(depCouronne_sf, crs = 2154)
paris_sf <- sf::st_transform(paris_sf, crs = 2154)

Pour être certain que les données et le territoire soient dans le même système de projection, il est possible de transformer ce dernier à l’aide de la fonction st_transform.

Territoires englobants, sélection des données à lisser

Pour éviter les effets de bord ==> sélectionner des données au-delà de notre zone d’intérêt (ici Paris intramuros)

Autrement dit, si on ne sélectionne que les ventes immobilières situées dans Paris intramuros, et qu’on lisse ce nuage de points, les zones situées à proximité du périphérique seront artificiellement peu denses en ventes immobilières. En effet, elles sont situées à proximité de zones artificiellement vides en ventes immobilières dans notre nuage de points, à savoir dans les communes limitrophes. Pourtant, dans la réalité, il y a bien des ventes réalisées à Malakoff, Vanves, Montrouge, etc. Et ces ventes influencent la densité lissée aux niveau de la Porte de Vanves.

Première méthode (sélection non-géométrique des ventes)

• filtrer les xy compris dans un grand rectangle englobant Paris intramuros
• Très efficace computationnellement
• Mais requiert de construire au préalable un rectangle adapté…

Ici, nous allons sélectionner les ventes immobilières appartenant à un grand rectangle englobant Paris intramuros.

L’idée est de pouvoir réduire la quantité de ventes à considérer avant le lissage, uniquement en filtrant les coordonnées x et y comprises dans le rectangle. Cette méthode est très efficace d’un point de vue calculatoire, mais requiert de construire au préalable un rectangle adapté.

• Mise en forme de la couche “territoire” : sélection des variables et renommage
• Création d’une bbox autour du territoire

paris_sf$nom <- "territoire"
paris_sf <- paris_sf[,c("nom","geometry")]

bbox <- sf::st_bbox(paris_sf)
bbox

   xmin    ymin    xmax    ymax 
 643076 6857499  660897 6867034 

Création d’un buffer de la bbox, avec une marge de 2000 mètres

marge <- 2000
bufferBbox <- bbox
bufferBbox[["xmin"]] <- bufferBbox[["xmin"]]-marge
bufferBbox[["xmax"]] <- bufferBbox[["xmax"]]+marge
bufferBbox[["ymin"]] <- bufferBbox[["ymin"]]-marge
bufferBbox[["ymax"]] <- bufferBbox[["ymax"]]+marge

Schéma explicatif des zones sélectionnées :

code du schema

# Petit rectangle vectoriel
bbox_sf = st_sf(nom="bbox", geometry = st_as_sfc(bbox), crs = 2154)

# Grand rectangle vectoriel
bufferBbox_sf = st_sf(nom="buffer_bbox", geometry = st_as_sfc(bufferBbox), crs = 2154)

# Options globales pour les cartes avec mapview
mapviewOptions(
  basemaps = c(
    "OpenStreetMap",
    "CartoDB.Positron",
    "CartoDB.DarkMatter",
    "Esri.WorldImagery"
  )
)
# Cartographie pédagogique avec mapview
m <- mapview(paris_sf ,col.regions= "#26cce7")+
  mapview(bufferBbox_sf %>% st_cast("MULTILINESTRING"),color="#FFC300",lwd=6)+
  mapview(bbox_sf %>% st_cast("MULTILINESTRING"),color="#229954",lwd=6)

m

• en bleu, le territoire d’intérêt (Paris)
• en vert, la bbox du territoire d’intérêt (à savoir le plus petit rectangle englobant cette zone)
• en jaune, la bbox élargie avec une marge de 2km = zone tampon autour de la bbox permettant de prendre en compte des observations au-delà de la seule zone d’intérêt, et d’ainsi éviter les effets de bord au moment du lissage.

Filtrage des coordonnées sont comprises à l’intérieur du rectangle jaune.

• très efficace computationnellement (filtre de valeurs numériques)
• Pas de traitement vectorielle (plus chronophages).

Filtrer (numériquement) les logements dans le triangle jaune

# Repérer les coordonnées extrêmes, avec une marge (ici 2km)
bufferBbox

   xmin    ymin    xmax    ymax 
 641076 6855499  662897 6869034 

xMin <- bufferBbox["xmin"]
xMax <- bufferBbox["xmax"]
yMin <- bufferBbox["ymin"]
yMax <- bufferBbox["ymax"]

# Ne garder que les données dans le rectangle englobant, sans traitement vectoriel !
dfBase_filtre <- dfBase[dfBase$x >= xMin & dfBase$x <= xMax & dfBase$y >= yMin & dfBase$y <= yMax, ]
dim(dfBase_filtre)

[1] 24419     8

dim(dfBase)

[1] 34489     8

Seconde méthode : la sélection géométrique des ventes

• Utiliser nos données individuelles comme un ensemble de points géolocalisés
• et procéder à des intersections géographiques.

1. On transforme nos observations en points vectoriels ;

sfBase <- sf::st_as_sf(dfBase, coords = c("x", "y"), crs = 2154)

2. On crée une zone tampon (buffer) autour du territoire d’intérêt, avec une marge (ici 2 000m), sous la forme d’un objet sf vectoriel ;

buffer_sf <- st_buffer(paris_sf, dist = 2000)

Remarque: Pour la zone tampon, prendre une marge légèrement plus grande que le rayon de lissage envisagé.

3. On repère les observations comprises dans cette zone tampon par intersection géographique.

sfBase_filtre <- st_join(sfBase,buffer_sf,left=F)
sfBase_filtre

Simple feature collection with 22221 features and 7 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 641099.2 ymin: 6855538 xmax: 662814.2 ymax: 6868990
Projected CRS: RGF93 / Lambert-93
First 10 features:
   id_mutation date_mutation  type_local nombre_pieces_principales
3  2021-447025    2021-01-08 Appartement                         3
5  2021-447029    2021-01-05 Appartement                         2
6  2021-447030    2021-01-15 Appartement                         3
9  2021-447034    2021-01-18 Appartement                         4
10 2021-447035    2021-01-18 Appartement                         4
11 2021-447037    2021-01-07 Appartement                         4
14 2021-447040    2021-01-07 Appartement                         4
15 2021-447041    2021-01-04 Appartement                         4
16 2021-447042    2021-01-15 Appartement                         3
19 2021-447046    2021-01-14 Appartement                         1
   valeur_fonciere surface_reelle_bati        nom                 geometry
3           384000                  41 territoire POINT (648001.8 6866153)
5           407200                  24 territoire POINT (646929.9 6864730)
6          1040000                  90 territoire POINT (645132.7 6864646)
9           625000                  68 territoire POINT (647494.6 6866019)
10          721250                  94 territoire POINT (643891.9 6864587)
11         1035000                  94 territoire   POINT (648239 6866235)
14          652000                  62 territoire POINT (647837.9 6865935)
15          830000                  94 territoire POINT (643093.5 6864422)
16          407070                  55 territoire POINT (647482.2 6867646)
19          235000                  22 territoire POINT (647066.8 6866298)

• Potentiellement lourd d’un point de vue calculatoire.
• Avec des données volumineuses, faire au minimum un premier filtrage non géométrique.

Ci-dessous le code permettant d’avoir un aperçu de la zone, du buffer et de 2000 points tirés aléatoirement dedans.

code du schema

# Mise en forme de la couche buffer
buffer_sf$nom <- "buffer"

# Échantillon de 2000 observations dans le buffer
sfBase_sample <- sfBase_filtre[sample(1:nrow(sfBase_filtre),2000) ,]

# Cartographie pédagogique avec mapview
mapview(paris_sf ,col.regions= "#26cce7")+
  mapview(buffer_sf %>% st_cast("MULTILINESTRING"),color="#FFC300",lwd=6)+
  mapview(sfBase_sample,#col.regions = "black",alpha.regions=0.5,
          alpha=0,cex=2)

Carroyage de données

Avant de lisser les données ponctuelles, on peut souhaiter représenter ces données sous forme carroyée afin de se les approprier. Il faut pour cela :

1. Associer chaque point (= vente géolocalisée) au centroïde du carreau auquel il appartient (le territoire est découpé en carreaux de 200 mètres à partir de l’origine du référentiel)

iCellSize = 200 # Carreaux de 200m
points_carroyage <- dfBase_filtre # On repart de la base filtrée selon la première méthode
points_carroyage$x_centroide = points_carroyage$x -
  (points_carroyage$x %% iCellSize) + (iCellSize / 2)
points_carroyage$y_centroide = points_carroyage$y -
  (points_carroyage$y %% iCellSize) + (iCellSize / 2)

head(points_carroyage)

  id_mutation date_mutation  type_local nombre_pieces_principales
1 2021-447023    2021-01-08 Appartement                         3
2 2021-447024    2021-01-05 Appartement                         2
3 2021-447025    2021-01-08 Appartement                         3
4 2021-447027    2021-01-08 Appartement                         3
5 2021-447029    2021-01-05 Appartement                         2
6 2021-447030    2021-01-15 Appartement                         3
  valeur_fonciere surface_reelle_bati        x       y x_centroide y_centroide
1          480000                  64 647357.3 6868635      647300     6868700
2          345000                  43 644483.5 6867695      644500     6867700
3          384000                  41 648001.8 6866153      648100     6866100
4          261900                  44 647414.8 6868937      647500     6868900
5          407200                  24 646929.9 6864730      646900     6864700
6         1040000                  90 645132.7 6864646      645100     6864700

2. Agréger les données sur chaque centroïde de la grille. En d’autres termes, compter le nombre de ventes par carreau

points_centroides <- points_carroyage %>%
  group_by(x=x_centroide,y=y_centroide) %>% 
  count(name = "nbVentes")
head(points_centroides,3)

# A tibble: 3 × 3
# Groups:   x, y [3]
       x       y nbVentes
   <dbl>   <dbl>    <int>
1 641100 6857700        1
2 641100 6857900        5
3 641100 6858300        1

Remarque : Ces deux premières étapes ne sont pas intégrées dans des fonctions du package btb. Il faut donc s’inspirer du code ci-dessus pour reproduire le même type de carroyage. Il est envisagé qu’une future version du package btb simplifie le code permettant de réaliser un carroyage.

3. Passer d’une table de centroïdes à une table de carreaux vectoriels

Utilisation de dfToGrid du package btb, avec comme paramètres obligatoires.

• df : un tableau avec les colonnes x et y représentant les coordonnées des centroïdes de la grille ;
• sEPSG : une chaîne de caractères indiquant le code epsg du système de projection utilisé ;
• iCellSize : la taille des carreaux (longueur du côté, en mètres).

carreaux <- btb::dfToGrid(df = points_centroides, sEPSG = "2154", iCellSize = iCellSize)

4. Se restreindre au champ des carreaux intersectant Paris

carreaux <- carreaux %>% st_join(paris_sf,left=F)

On obtient le carroyage des ventes dans Paris intramuros

code du production de la carte

contourParis <- st_cast(paris_sf[,c("geometry")],"MULTILINESTRING")

mf_init(x=carreaux,theme = "agolalight")
mf_map(x = carreaux,
       type = "choro",
       var="nbVentes",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       leg_val_rnd = 1,
       add = TRUE)
mf_map(x = contourParis,
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Carroyage du nombre de ventes",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Remarque :

Le carroyage pourrait aussi être utilisé pour simplifier les données avant lissage si le calcul de lissage est trop long en raison d’un très grand nombre d’observations.

Lissage

Calcul de densité

• Lissage le plus simple à réaliser
• Quantité par unité de surface (un carreau)
• Ici : densité de ventes de logements dans la ville de Paris au cours de l’année 2021

Premier lissage de densité

Toujours pour éviter les effets de bord :

• Lissage sur une zone plus large que la zone d’intérêt
• Utilisation de dfBase_filtre (construite précédemment).

Premier lissage de densité

• Création une variable nbObsLisse qui vaudra 1 pour chaque observation
• et on lisse nbObsLisse avec la fonction kernelSmoothing du package btb.

dfBase_filtre$nbObsLisse <- 1
head(dfBase_filtre)

  id_mutation date_mutation  type_local nombre_pieces_principales
1 2021-447023    2021-01-08 Appartement                         3
2 2021-447024    2021-01-05 Appartement                         2
3 2021-447025    2021-01-08 Appartement                         3
4 2021-447027    2021-01-08 Appartement                         3
5 2021-447029    2021-01-05 Appartement                         2
6 2021-447030    2021-01-15 Appartement                         3
  valeur_fonciere surface_reelle_bati        x       y nbObsLisse
1          480000                  64 647357.3 6868635          1
2          345000                  43 644483.5 6867695          1
3          384000                  41 648001.8 6866153          1
4          261900                  44 647414.8 6868937          1
5          407200                  24 646929.9 6864730          1
6         1040000                  90 645132.7 6864646          1

Premier lissage de densité

Paramètres de kernelSmoothing :

• dfObservations : le tableau des données à lisser. Il doit nécessairement contenir une colonne x, une colonne y, et 1 à n colonnes numériques (variables à lisser) ;
• sEPSG : chaine de caractères indiquant le code epsg du système de projection utilisé ;
• iCellSize : un entier indiquant la taille des carreaux ;
• iBandwidth : un entier indiquant le rayon de lissage.

Premier lissage de densité

Lissage et enregistrement dans sfCarrLiss

dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
sfCarrLiss <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 400)

Premier lissage de densité

Visualisation du résultat :

head(sfCarrLiss)

Simple feature collection with 6 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 646000 ymin: 6855400 xmax: 647800 ymax: 6855600
Projected CRS: RGF93 / Lambert-93
       x       y nbObsLisse                       geometry
1 646100 6855500   1.712851 POLYGON ((646000 6855600, 6...
2 646500 6855500   1.261316 POLYGON ((646400 6855600, 6...
3 647100 6855500   2.301557 POLYGON ((647000 6855600, 6...
4 647300 6855500   2.808277 POLYGON ((647200 6855600, 6...
5 647500 6855500   2.395913 POLYGON ((647400 6855600, 6...
6 647700 6855500   1.949819 POLYGON ((647600 6855600, 6...

Premier lissage de densité

Nombre de lignes lissées

nrow(sfCarrLiss)

[1] 4107

Premier lissage de densité

code

mf_init(x=paris_sf,theme = "agolalight")
mf_map(x = st_cast(sfCarrLiss[,c("geometry")],"LINESTRING"), 
       lwd=1,add=T)
mf_map(x = contourParis, 
       lwd=8,
       col="wheat",add = TRUE)
mf_layout(title = "KernelSmoothing génère une grille carroyée",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Premier lissage de densité

Remarque :

• la grille lissée épouse le périmètre des points en entrée.
• Elle va même un peu au-delà (voir infra).

Premier lissage de densité

Attention, la fonction retourne une erreur en cas de :

• présence d’une variable non-numérique ;
• valeur(s) absente(s) dans les colonnes x ou y.

Premier lissage de densité

La fonction de lissage classique du package btb est conservative.

# Nombre de ventes dans la petite couronne
sum(dfBase_filtre$nbObsLisse)

[1] 24419

# Après lissage, nombre lissé de ventes dans les carreaux
sum(sfCarrLiss$nbObsLisse) 

[1] 24419

Premier lissage de densité

Remarque sur le secret :

• Toujours vérifier son respect avant publication de cartes !
• Par exemple, avec données fiscales, pas de carreaux comportant moins de 11 observations, même lissées…

sfCarrLiss <- sfCarrLiss[sfCarrLiss$nbObsLisse >= 11, ]

Premier lissage de densité (représentation)

On cartographie les carreaux :

• On réalise une carte de type choroplèthe ;
• On discrétise de la densité lissée (ici : quantiles) ;
• Avec 5 classes (c’est-à-dire 5 couleurs).

Premier lissage de densité (représentation)

code

# Filtrage des carreaux lissés intersectant la ville de Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf[,"geometry"],left=F)

# Carte lissée
mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 400m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Remarque : la variable lissée nbObsLiss correspond au nombre de ventes (observations) par carreau (ici de 200 mètres). Ainsi, il est possible d’obtenir une densité au km² en multipliant la variable nbObsLiss par 25 dans le cas présent.

Premier lissage de densité (recapitulatif)

On souhaite obtenir le nombre de ventes lissé, en 2021 dans Paris intramuros :

• On sélectionne les ventes au-delà de Paris intramuros (effets de bord)
• On lisse grâce à kernelSmoothing
• On obtient une grille carroyée bien plus vaste que Paris intramuros
• On filtre uniquement les carreaux dans Paris
• On cartographie la variable lissée (nbObsLissee)

Faire varier le rayon de lissage

Taille optimale du rayon de lissage ?

Inconnu a priori…

Compromis entre précision et généralisation

On teste avec :

• 600m
• puis 1km

Faire varier le rayon de lissage (600m)

code lissage

dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
sfCarrLiss <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 600)
# Filtrage des carreaux lissés dans Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf[,"geometry"],left=F)

code carto

# Carte lissée
mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 600m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Faire varier le rayon de lissage (1000m)

code lissage

dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
sfCarrLiss <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 1000)
# Filtrage des carreaux lissés dans Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf[,"geometry"],left=F)

code carto

# Carte lissée
mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       lwd=1,
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 1000m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Faire varier le rayon de lissage

• Rayon grand : aspects structurels des données
• Rayon petit : spécificités locales des données.

Le choix revient donc au statisticien-géographe, en fonction de sa connaissance des données et des objectifs recherchés.

Sans la grille apparente

code carto

mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       border = NA, # C'est ici que ça se passe
       # lwd=1,
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage avec rayon de 1000m, sans grille",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Eviter l’effet granuleux

Arètes des carreaux visibles

Idée : réduite la taille des carreaux (=> 50m)

Attention au temps de calcul !

(et on en profite aussi pour faire varier le rayon de lissage => 1km)

Eviter l’effet granuleux

code lissage

dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
sfCarrLiss <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 50, 
                                    iBandwidth = 1000)
# Filtrage des carreaux lissés dans Paris
sfCarrLiss_paris <- sfCarrLiss %>% st_join(paris_sf[,"geometry"],left=F)

code carto

mf_init(x=sfCarrLiss_paris,theme = "agolalight")
mf_map(x = sfCarrLiss_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = "quantile",
       nbreaks = 5,
       border = NA,
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage avec pas de 50m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Illustration des effets de bord

Que se passe-t-il si on lisse uniquement les ventes réalisées dans Paris intramuros ?

Comparons avec les mêmes paramètres :

• Rayon de lissage de 2 000 mètres
• Taille des carreaux de 200m
• Bornes de discrétisation communes (comparabilité)

Sans effets de bord

code lissage

dataLissage <- dfBase_filtre[,c("nbObsLisse","x","y")]
sfCarLiss <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 2000)
sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf,left = F)

code carte

mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = c(0,3,5,7,9,17),
       # nbreaks = 5,
       border = NA,
       leg_val_rnd = 1,
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Sans effets de bord, avec R=2000m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Avec effets de bord

code lissage

sfBase_intramuros <- sfBase %>% st_join(paris_sf,left=F)
dfBase_intramuros <- sfBase_intramuros %>% 
  mutate(x=st_coordinates(geometry)[,1],
         y=st_coordinates(geometry)[,2],
         nbObsLisse=1) %>% 
  st_drop_geometry() %>% 
  select(nbObsLisse,x,y)

sfCarLiss_intramuros <- btb::kernelSmoothing(dfObservations = dfBase_intramuros, 
                                    sEPSG = "2154",
                                    iCellSize = 200, 
                                    iBandwidth = 2000)

sfCarLiss_intramuros_paris <- sfCarLiss_intramuros %>%
  st_join(paris_sf,left = F)

code carte

mf_init(x=sfCarLiss_intramuros_paris,theme = "agolalight")
mf_map(x = sfCarLiss_intramuros_paris, 
       type = "choro",
       var="nbObsLisse",
       breaks = c(0,3,5,7,9,17),
       # nbreaks = 5,
       border = NA,
       leg_val_rnd = 1, 
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Avec effets de bord, avec R=2000m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Illustrons des effets de bord (conclusion)

Les effets de bord se matérialisent à la périphérie de Paris, où les densités lissées sont artificiellement plus faibles que quand on lisse en prenant une marge.

Les grilles automatiques de btb

Par défaut : la grille produite automatiquement par btb::kernelSmoothing dépasse les limites de la zone d’étude choisie.

Prenons l’exemple précédent du lissage sur les ventes à l’intérieur de Paris intramuros.

Les grilles automatiques de btb

code carte

mapview(sfCarLiss_intramuros %>%
          st_cast("MULTILINESTRING"), lwd = 2,
        layer.name="Grille automatique") +
  mapview(contourParis,color="red",lwd=10) 

La grille obtenue dépasse Paris intramuros et contient :

• l’ensemble des carreaux contenant au moins une vente…
• + 4 carreaux limitrophes (règle automatique - cf. formation « Comment utiliser les outils de l’analyse urbaine »)

Les grilles automatiques de btb

Modification de cette règle avec le paramètre iNeighbor.

Par exemple :

• iNeighbor=2 => 2 carreaux limitrophes aux carreaux contenant au moins une vente.
• iNeighbor=0 => grille uniquement constituée des carreaux contenant au moins une vente
  • attention à l’effet « gruyère »

Les grilles automatiques de btb

code lissage

sfCarLiss_intramuros <- btb::kernelSmoothing(
  dfObservations = dfBase_intramuros, 
  sEPSG = "2154",
  iCellSize = 200, 
  iBandwidth = 2000,
  iNeighbor = 0)

code carte

mapview(sfCarLiss_intramuros %>% st_cast("MULTILINESTRING"),lwd=2,layer.name="Grille automatique") +
  mapview(contourParis,color="red",lwd=10) 

Les grilles automatiques de btb

On peut souhaiter utiliser une grille de lissage particulière et adaptée à un territoire donné.

Par exemple, sur une ville cotière, on ne veut :

• ni de carreau dans la mer
• ni l’effet « gruyère » du paramétrage iNeighbor=0.

Dans ce cas :

• il convient de la confectionner soi-même
• et de la renseigner dans le paramètre dfCentroids de la fonction kernelSmoothing.

Calcul de moyenne

Dans cette section, nous allons nous intéresser au prix moyen des logements vendus à Paris en 2021.

Remarque :

Une moyenne n’est pas le lissage du rapport liss(variable / nbObsLisse) mais le rapport des lissages liss(variable) / liss(nbObs).

Calcul de moyenne : comment faire ?

Il faut :

1. Lisser les prix de ventes des biens vendus
2. Lisser le nombre de biens vendus (exactement comme précédemment)
3. Faire le ratio des deux précédents lissages pour chaque carreau de la grille produite.

Calcul de moyenne : comment faire ?

code lissage

dataLissage <- dfBase_filtre[,c("nbObsLisse","valeur_fonciere","x","y")] 

# Lissage
sfCarLiss <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 100, 
                                    iBandwidth = 800)
sfCarLiss

code carte

# Calculer le taux
sfCarLiss$prixMoyen <- sfCarLiss$valeur_fonciere / sfCarLiss$nbObsLisse

# Cartographie
sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf,left=F)
mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris, 
       type = "choro",
       var="prixMoyen",
       breaks = "quantile",
       nbreaks = 5,
       border=NA,
       leg_val_rnd = 1, 
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Lissage des prix de vente avec un rayon de 800m",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Calcul de moyenne : resultats

Les prix moyens sont particulièrement élevés dans le centre et l’Ouest de la capitale

Remarque :

le lissage en moyenne est potentiellement sensible aux valeurs atypiques (contrairement au lissage quantile)

Calcul de taux

Taux de grands logements

Nous nous intéressons ici à la proportion de ventes portant sur des grands logements (comportant plus de 4 pièces principales).

Création de la variable quatrePieces (indicatrice)

dfBase_filtre$quatrePieces <- ifelse(
  dfBase_filtre$nombre_pieces_principales >= 4, 1, 0)

Taux de grands logements

• Lissage des variables quatrePieces et nbObsLisse
• Ratio des deux variables lissées

# Lissage
dataLissage <- dfBase_filtre[,c("nbObsLisse","quatrePieces","x","y")]
sfCarLiss <- btb::kernelSmoothing(dfObservations = dataLissage, 
                                    sEPSG = "2154",
                                    iCellSize = 100, 
                                    iBandwidth = 800)
# Calculer le taux
sfCarLiss$txQuatrePieces <- sfCarLiss$quatrePieces / sfCarLiss$nbObsLisse

# Intersection géographique avec Paris
sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf,left=F)

Taux de grands logements

Code

mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris, 
       type = "choro",
       var="txQuatrePieces",
       breaks = "quantile",
       nbreaks = 5,
       border=NA,
       leg_val_rnd = 1, 
       add = TRUE)
mf_map(x = contourParis, 
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Logements avec 4 pièces ou plus (rayon de 800m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Taux de grands logements

Fort logiquement, la carte des prix moyens et la carte des grands logements ont des similitudes importantes.

On peut alors avoir envie de lisser le prix au mètre-carré.

Prix au m²

Il convient de lisser séparément les prix de vente (numérateur) et le nombre de mètres-carrés (dénominateur) des logements vendus.

Autrement, si on lisse le prix au m² de chaque logement vendu, on sur-pondère artificiellement les prix au m² des petits logements (car l’unité statistique devient alors le logement, et non le m² vendu).

Prix au m²

code lissage

dataLissage <- dfBase_filtre[,c("surface_reelle_bati","valeur_fonciere","x","y")]

sfCarLiss <- btb::kernelSmoothing(dfObservations = dataLissage,
                                    sEPSG = "2154",
                                    iCellSize = 100,
                                    iBandwidth = 800)
# Calculer le prix au m²
sfCarLiss$prixM2 <- sfCarLiss$valeur_fonciere / sfCarLiss$surface_reelle_bati

code carte

sfCarLiss_paris <- sfCarLiss[unlist(st_intersects(paris_sf,sfCarLiss)),]

mf_init(x = sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
       type = "choro",
       var="prixM2",
       breaks = "quantile",
       nbreaks = 5,
       border=NA,
       leg_val_rnd = 1, 
       add = TRUE)
mf_map(x = st_cast(paris_sf[,c("geometry")],"MULTILINESTRING"),
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Prix au m² (rayon de 800m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Calcul de quantiles géographiquement pondérés

Lissage quantile : pourquoi ?

• Moins sensible aux valeurs extrêmes ;
• Obtenir des indicateurs de dispersion lissés ;
  • ex : écart interquantile

Lissage quantile : comment ?

• Toujours avec la fonction kernelSmoothing
• En ajoutant un paramètre, vQuantiles, qui contient la liste des quantiles souhaités

Lissage quantile : rappel

Contrairement au lissage classique :

• il n’est pas conservatif
• la fonction crée :
  • une colonne nbObs indiquant le nombre réel (non lissé) d’observations ayant contribué au calcul du carreau.
  • pour chaque variable, autant de colonnes qu’il y a de quantiles souhaités

Lissage quantile : exemple

Nous allons calculer le 1er décile, la médiane et le 9ème décile du prix de vente des logements à Paris.

dataLissage <- dfBase_filtre[,c("valeur_fonciere","x","y")]
sfCarLiss <- btb::kernelSmoothing(dfObservations = dataLissage,
                                    sEPSG = "2154",
                                    iCellSize = 100,
                                    iBandwidth = 1500,
                                    vQuantiles = c(0.1, 0.5, 0.9))

Lissage quantile : exemple

Visualisons les première lignes de la grille lissée

head(sfCarLiss)

Simple feature collection with 6 features and 6 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 651500 ymin: 6854700 xmax: 652100 ymax: 6854800
Projected CRS: RGF93 / Lambert-93
  nbObs valeur_fonciere_01 valeur_fonciere_05 valeur_fonciere_09      x       y
1    57             182000             310000             610000 651550 6854750
2    53             182000             310000             620000 651650 6854750
3    52             182000             310000             620000 651750 6854750
4    55             182000             321000             620000 651850 6854750
5    48             182000             324000             620000 651950 6854750
6    55             182000             324000             620000 652050 6854750
                        geometry
1 POLYGON ((651500 6854800, 6...
2 POLYGON ((651600 6854800, 6...
3 POLYGON ((651700 6854800, 6...
4 POLYGON ((651800 6854800, 6...
5 POLYGON ((651900 6854800, 6...
6 POLYGON ((652000 6854800, 6...

Lissage quantile : exemple

Cartographie du 1er décile de prix des ventes

code carte

sfCarLiss_paris <- sfCarLiss %>% st_join(paris_sf,left = F)
mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
       type = "choro",
       var="valeur_fonciere_01",
       breaks = "quantile",
       nbreaks = 5,
       border=NA,
       leg_val_rnd = 0,
       add = TRUE)
mf_map(x = contourParis,
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Premier décile des prix de vente (rayon de 1500m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Lissage quantile : exemple

Cartographie de la médiane des prix des ventes

code carte

mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
       type = "choro",
       var="valeur_fonciere_05",
       breaks = "quantile",
       nbreaks = 5,
       border=NA,
       leg_val_rnd = 0,
       add = TRUE)
mf_map(x = contourParis,
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Médiane des prix de vente (rayon de 1500m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

Lissage quantile : exemple

Cartographie du rapport interdécile des prix des ventes

code carte

sfCarLiss_paris$rapp_interdecile <- sfCarLiss_paris$valeur_fonciere_09 / sfCarLiss_paris$valeur_fonciere_01

# Cartographie
mf_init(x=sfCarLiss_paris,theme = "agolalight")
mf_map(x = sfCarLiss_paris,
       type = "choro",
       var="rapp_interdecile",
       breaks = "quantile",
       nbreaks = 5,
       border=NA,
       leg_val_rnd = 0,
       add = TRUE)
mf_map(x = contourParis,
       lwd=4,
       col="black",add = TRUE)
mf_layout(title = "Rapport interdéciles des prix de vente (rayon de 1500m)",
          credits = "Insee-DSAU, DGFiP, Etalab, IGN, mapsf")

(attention, il est ici question des prix de vente, et non des prix au m²).

Lissage quantile : remarque importante

Enfin, pour que le quantile lissé ait du sens, il faut qu’un nombre conséquent d’observations ait participé à sa construction. Ainsi, il est conseillé d’utiliser un rayon de lissage suffisamment élevé.

summary(sfCarLiss_paris$nbObs)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     13     725    1040    1106    1440    2917 

sfCarLiss_paris %>% st_drop_geometry() %>% group_by(nbObs<100) %>% count()

# A tibble: 2 × 2
# Groups:   nbObs < 100 [2]
  `nbObs < 100`     n
  <lgl>         <int>
1 FALSE         10200
2 TRUE            229

Indicateurs sur zonage à façon

Objectif : agréger des données sur un zonage particulier

2 cas :

1. données ponctuelles
2. données carroyées

Données ponctuelles

Calculons le nombre et le prix moyen des transactions immobilières en 2021 dans le triangle d’or.

Pour information, ce triangle a été créé manuellement sur le site du Géoportail.

Données ponctuelles

1. Import et visualisation de la zone à façon « Triangle d’or ». Ce fichier géographique est au format .kml et se lit sans problème avec la fonction sf::st_read.

chemin_file <- paste0(url_bucket,bucket,"/r-lissage-spatial/triangle_or.kml")
triangleOr <- st_read(chemin_file)

Données ponctuelles

On visualise le triangle :

mapview(triangleOr,col.regions="#ffff00")

Données ponctuelles

2. Intersection géographique entre des points (les transactions) et un polygone (le quartier). Attention au système de projection.

st_crs(triangleOr)[["input"]] # Système de projection de triangleOr

[1] "WGS 84"

triangleOr <- st_transform(triangleOr,crs = 2154) # Transformation en Lambert93

# Intersection géographique avec les transactions
transac_triangle <- sfBase %>% st_join(triangleOr[,"geometry"],left=F)

head(transac_triangle,4) # Visualisation de points retenus

Simple feature collection with 4 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 648814.5 ymin: 6863055 xmax: 649194 ymax: 6863601
Projected CRS: RGF93 / Lambert-93
      id_mutation date_mutation  type_local nombre_pieces_principales
18629 2021-489256    2021-01-15 Appartement                         5
18676 2021-489336    2021-01-25 Appartement                         3
18677 2021-489337    2021-01-28 Appartement                         4
18698 2021-489369    2021-02-05 Appartement                         6
      valeur_fonciere surface_reelle_bati                 geometry
18629         3300000                 180 POINT (648918.9 6863145)
18676          590000                  53 POINT (648903.5 6863601)
18677         1511345                  65   POINT (649194 6863377)
18698         2450000                 133 POINT (648814.5 6863055)

Données ponctuelles

Cartographie des points dans le triangle

mapview(transac_triangle)+mapview(triangleOr,col.regions="#ffff00")

Données ponctuelles

3. Calculer du nombre et du prix moyen des transactions immobilières réalisées en 2021 dans ce quartier

# Nombre de transactions en 2021 : 
nrow(transac_triangle)

[1] 20

# Prix moyen
mean(transac_triangle$valeur_fonciere)

[1] 1637256

Remarque : On ne voit que 12 points sur la cartes, alors qu’on a trouvé 20 transactions dans le triangle : ceci s’explique car certaines transactions sont géolocalisées au même endroit (appartements situés dans un même immeuble par exemple).

Données carroyées

Objectif : obtenir la proportion de logements sociaux dans un rayon de 1km autour de la Direction générale de l’Insee.

On utilise les données carroyées à 200 mètres disponibles sur insee.fr.

Remarque : un programme généralisant cet exercice est disponible sur insee.fr, en accompagnement des données carroyées.

Données carroyées

1. Chargement des données carroyées de Filosofi 2015 à 200m, uniquement en région parisienne, en ne conservant que 6 variables dont le nombre de logements sociaux.

# Fonction permettant de faire le chargement souhaité
st_read_maison <- function(chemin_tab){
  requete <- "SELECT IdINSPIRE,Depcom,I_est_cr,Men, Log_soc, geom
            FROM Filosofi2015_carreaux_200m_metropole
            WHERE SUBSTR(Depcom, 1, 2) IN ('75','92','93','94') "
  sf::st_read(chemin_tab, query = requete)
}

# Chargement des données
chemin_file <- paste0(url_bucket, bucket, 
                    "/r-lissage-spatial/Filosofi2015_carreaux_200m_metropole.gpkg")
carreaux <-  st_read_maison(chemin_file)
head(carreaux)

Données carroyées

head(carreaux)

Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 652623.2 ymin: 6858353 xmax: 655268.9 ymax: 6866380
Projected CRS: RGF93 / Lambert-93
                       IdINSPIRE Depcom I_est_cr Men Log_soc
1 CRS3035RES200mN2893400E3763200  75119        0 990     937
2 CRS3035RES200mN2890000E3762200  75111        0 926      80
3 CRS3035RES200mN2893400E3762400  75119        0 508     323
4 CRS3035RES200mN2891800E3763600  75119        0 633       0
5 CRS3035RES200mN2890800E3763200  75120        0 349     188
6 CRS3035RES200mN2885800E3760600  75113        0 751     344
                            geom
1 MULTIPOLYGON (((654521.9 68...
2 MULTIPOLYGON (((653843.4 68...
3 MULTIPOLYGON (((653725.1 68...
4 MULTIPOLYGON (((655069.7 68...
5 MULTIPOLYGON (((654764.7 68...
6 MULTIPOLYGON (((652641.8 68...

Données carroyées

2. Créer un disque de 1km autour du White

white <- st_sfc(st_point(c(649218.36,6857569.14)),crs=2154)
buffer_white <- st_buffer(white,1000) %>% st_as_sf()

Données carroyées

3. Déterminer, pour cette zone, l’ensemble des carreaux qui la recouvrent

Remarque : L’ensemble de carreaux étant en général plus large que la zone, les agrégats obtenus seront des estimations des valeurs réelles, plus ou moins précises. Le problème ne se pose pas si vous travaillez directement sur des données disponibles au niveau « individu statistique » (avec des coordonnées x,y).

carreaux_select <- carreaux %>% st_join(buffer_white,left=F)

Données carroyées

À noter que ces carreaux sont « penchés » : ceci est la résultante de leur reprojection en Lambert 93. En effet, ces carreaux ont été produits en projection LAEA (standard européen), ce qui leur donne cet aspect une fois reprojetés dans le standard français.

mapview(carreaux_select)+
  mapview(buffer_white,color="#ffff00",lwd=10,alpha.regions=0,legend=F)

Données carroyées

4. Calculer l’indicateur souhaité sur la zone par somme des données par carreau.

tx_logSoc_white <- sum(carreaux_select$Log_soc)/sum(carreaux_select$Men)
cat("On compte ",
    round(tx_logSoc_white*100),
    "% de logements sociaux dans un rayon de 1km autour du white")

On compte  33 % de logements sociaux dans un rayon de 1km autour du white

Conclusion

Vous êtes désormais capables de carroyer et lisser toutes les données que vous aurez à votre disposition !

Autre logiciels

• Python : libraryBTBpy
• Qgis grâce à un plug-in développé par Lionel Cacheux (DR Insee Grand Est). Tutoriel à destinations des agents de l’Insee disponible sur l’interface AUS (V:/PSAR-AU/Formation Comment utiliser les outils AU/2019/Séquence 6 - Outils/Formation_AU_Sequence_6_TP_partie_1.odt)
• ALICE, application intranet Insee développée par le PSAR Analyse urbaine (non maintenue à ce jour)

Références :

• Manuel d’analyse spatiale de l’Insee
• Introduction à la géomatique pour le statisticien, Document de travail, Insee
• Github du package BTB